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Abstract 

Background: In photoacoustic imaging (PAI), the reduction of scanning time is a 
major concern for PAI in practice. A popular strategy is to reconstruct the image from 
the sparse-view sampling data. However, the insufficient data leads to reconstruction 
quality deteriorating. Therefore, it is very important to enhance the quality of the 
sparse-view reconstructed images. 

Method: In this paper, we proposed a joint total variation and /.p-norm (hZ-Lp) based 
image reconstruction algorithm for PAI. In this algorithm, the reconstructed image is 
updated by calculating its total variation value and i.p-norm value. Along with the 
iteration, an operator-splitting framework is utilized to reduce the computational cost 
and the Barzilai-Borwein step size selection method is adopted to obtain the faster 
convergence. 

Results and conclusion: Through the numerical simulation, the proposed algorithm is 
validated and compared with other widely used PAI reconstruction algorithms. It is 
revealed in the simulation result that the proposed algorithm may be more accurate 
than the other algorithms. IVloreover, the computational cost, the convergence, the 
robustness to noises and the tunable parameters of the algorithm are all discussed 
respectively. We also implement the JM-Lp algorithm in the in-vitro experiments to 
verify its performance in practice. Through the numerical simulations and in-vitro 
experiments, it is demonstrated that the proposed algorithm enhances the quality 
of the reconstructed images with faster calculation speed and convergence. 

Keywords: Photoacoustic imaging. Image reconstruction, Total variation, /.p-norm, 
Nonconvex optimization 



Introduction 

Photoacoustic imaging (PAI), also known as optoacoustic tomography (OAT) or ther- 
moacoustic tomography (TAT), is a novel hybrid biomedical imaging modality which 
combines the strengths of both optical and ultrasound imaging [1-6]. Due to its non- 
ionizing nature, it has been considered as a promising imaging technique and developed 
rapidly during the past decade. PAI reveals physiologically specific optical absorption 
contrast of the biological tissues, which has great potential in clinic applications such as 
early tumor detection [7,8], vessel imaging [9,10] and brain imaging [11]. 

PAI is developed based on the photoacoustic effect [1,2], which is a process describ- 
ing that the imaging tissues absorb the laser energy and convert it into acoustic waves. 
In this paper, we focus on the computed-tomographic PAI. In this kind of imaging 
mode, a laser pulse is used to illuminate the imaging tissues from the top. Some of the 
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laser energy is absorbed and converted into heat, leading to thermoelastic expansion 
and thus wideband ultrasonic wave emission. The generated photoacoustic signals are 
then detected by a scanning ultrasound transducer or a transducer array to form im- 
ages. Based on these detected signals, the optical absorption deposition of the imaging 
tissues can be calculated by using an image reconstruction algorithm. 

In PAI, the reconstruction algorithms have become the vital factor of imaging quality. 
The PAI reconstruction result benefits a lot from a stable, accurate and efficient algo- 
rithm. A variety of analytical image reconstruction algorithms have been developed. Re- 
construction algorithms based on the inverse spherical radon transform have been 
proposed in both the time-domain [12,13] and the frequency-domain [14,15]. The fil- 
tered back-projection (FBP) algorithm proposed by Xu et al. is the most popular algo- 
rithm due to its accuracy and convenience [16,17]. The deconvolution reconstruction 
algorithm proposed by Zhang et al. has specific advantage under the circumstance of 
limited-angle sampling and heterogeneous acoustic medium [18,19]. Several investiga- 
tions have been made to propose the algorithms in plane geometries for imaging with 
the linear array of transducer [20,21]. The analytical reconstruction methods have 
advantage in the computational cost and implementation convenience. However, the 
analytical algorithms fail to keep effective when the sampling points are sparse. The 
ignorance of measurements noises leads to the severely quality decline in the noisy 
situation. Those drawbacks above limit the applications of the analytical algorithms and 
impair their performance. Then the iterative image reconstruction methods are pro- 
posed to overcome these shortcomings and enhance the image quality of PAI. 

The iterative image reconstruction methods usually build up a model to describe the 
relationship between the detected photoacoustic signals and the optical absorption de- 
position. So they are also called the model-based algorithms. Most of them calculate 
the optical absorption deposition iteratively to get the final reconstructed image. With 
proper optimization condition setup, the model-based methods can provide a more ac- 
curate and robust image reconstruction compared to the analytical ones [22,23]. Many 
methods that proved to be useful in other aspects have been adopted in PAI recon- 
struction as an optimization condition of the model-based methods. Some algorithms 
focus on the compensation of the acoustic inhomogeneous phenomenon [24,25]. Jose 
et al. proposes an iterative approach that takes the speed-of-sound of subject into ac- 
count. They acquire the 2D speed-of-sound distributions and use this speed-of-sound 
map in their reconstruction algorithm [24]. Huang et al. develop and establish a full- 
wave iterative reconstruction approach in PAI to deal with the acoustic inhomogeneous 
and acoustic attenuation problem [25]. The compressed sensing has been involved in 
PAI reconstruction aiming to reduce the measurements and accelerate the data acquisi- 
tion [26,27]. The model-based algorithm proposed by Rosenthal et al. recovers the 
image in the wavelet domain with a different strategy [28]. Meng et al. develop a com- 
pressed sensing framework by using partially known support [29,30]. The reported re- 
sults show some improvement of image qualities. The total variation (TV) coefficient is 
always used to de-noise the image. Some algorithm are proposed by using the total 
variation minization to PAI image reconstruction [31-34]. Yao et al. propose the total 
variation minization (TVM) with the TV coefficient involved in the finite elment 
method to enhance the image quality and overcome the limit-angle problem [31,32]. 
An adaptive steepest-descent-projection onto convex sets (ASD-POCS) is proposed by 
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Wang et al. with the TV utilized in the iteration [33]. They investigate and employ the 
TV-based iterative image reconstruction algorithms in three-dimensional PAL Zhang 
et al. utilizes the TV coefficient along with the gradient descent method in PAI recon- 
struction to propose the Total variation based gradient descent (TV-GD) algorithm 
[34]. The TV-GD method is reported to be stable and efficient under the sparse-view 
circumstance for PAI reconstruction. From the discussion noted above, it can be de- 
duced that the iterative algorithms for PAI reconstruction have advantage in recon- 
struction qualities and robustness to noise. Now the reduction of scanning time is the 
main concern of PAI. A popular strategy is to reconstruct the image from sparse-view 
sampling data. Also there exists photoacoustic imaging system which can image the 
whole area with one laser exposure. These systems usually have large amount of trans- 
ducers around the imaging area. With the help of sparse-view photoacoustic imaging 
reconstruction method, the transducer amount can be reduced. This reduction benefits 
the system from two main aspects. First, it helps to manage the system complexity to a 
lower level. The lower complexity system is more stable and easier to maintain. Second, 
this reduction also means the reduction of data scale. The data scale reduction can 
make the acquisition process more simple and flexible. Besides these two aspects, it is 
also worth to mention that it reduces the cost of the whole system. All those men- 
tioned above is very important for further clinical applications. It is very important to 
develop a sparse-view imaging system. In this situation, the qualities of the iterative re- 
constructed images have room for improvement. Take the TV-GD method for example, 
it is reported to be an efficient and high-quality algorithm in sparse-view situation. But 
the paintinglike artifacts emerge and some detail information is lost in the extremely 
sparse-view reconstruction. 

The compressed sensing theories have been adopted in PAI reconstruction, in which 
the Li-norm of the signal is minimized to obtain the reconstructed image. Recently, 
Chartrand reported that by replacing the ii-norm with the i^-norm {Q<p< 1), accur- 
ate reconstruction is possible with substantially fewer measurements [35]. This noncon- 
vex optimization setting has been successfully applied to Magnetic Resonance Imaging 
(MRI) image reconstruction [36,37]. The results show that the algorithm with L^-norm 
can provide accurate reconstruction image with fewer measurements comparing to the 
Zi-norm based algorithms. To another dimension of this optimization problem, seve- 
ral algorithms have been proposed to get a better performance in image reconstruc- 
tion through jointly minimizes the TV value and Lj-norm value [38,39] in MRI image 
reconstruction. 

In this paper, we present a novel algorithm to the problem of reconstructing the image 
from sparse- view data in PAI. The algorithm is based on the jointly minimization of total 
variation and nonconvex L^-norm (TV-L^). The reconstructed image is updated by calcu- 
lating its joint total variation value and L^-norm value. The operator-splitting framework 
is used to reduce the computational cost, and the Barzilai-Borwein step size selection 
method is adopted to obtain the faster convergence. Through the numerical simulation, 
the image reconstruction in the case of insufficient sampling data was accomplished. The 
reconstruction result is compared with several other algorithms including the FBP [16], 
the Li-norm [27] and the TV-GD method [34]. The computational cost and the conver- 
gence of the proposed algorithm is also discussed and compared with other algorithms. 
The numerical simulations also cover the robustness to the noise and the tunable 
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parameter investigation. Through the numerical simulations and in-vitro experiments, it 
is demonstrated that the proposed algorithm enhances the quality of the reconstructed 
images with the faster calculation speed and convergence. It's worthwhile to mention that 
like Ref. [27] and other iteration method, we also used a projection matrix to connect the 
acoustic pressure measurements with the reconstructed image. But there are some imple- 
mentation differences between our method and that one. We both use an intermediate 
variable to simplify our equations. Ref [27] used the velocity potential as the intermediate 
variable and we used a linear integration of the initial pressure along an arc whose center 
is the position of the ultrasound sensor and with a certain radius ct. The Ref [27] used a 
sparsifying matrix and minimized the Li-norm in sparsifying domain to get the recon- 
struction. We used the information from sparsifying domain and piecewise continuous 
behavior to reconstruct the image. Also, we adapted the /i-norm minimization into the al- 
gorithm, so it can be a more accurate algorithm in sparse-view PAL 

The main contribution of this paper is to develop a novel algorithm for solving the 
problem of reconstructing the image from sparse-view data in PAI. Our contributions 
are threefold. First, we include the nonconvex optimization into the PAI reconstruction. 
This nonconvex optimization setting can provide more stable and accurate result under 
the sparse-view situation. Second, we combine the nonconvex optimization with TV 
minimization. The combined method is able to reconstruct more detailed image with 
sharp edge. Finally, we implement the Barzilai-Borwein method accelerates the recon- 
struction speed and improves the convergence considerably. 

This paper is organized as follows. 'Theory and method' describes the theory of the 
proposed algorithm. The numerical simulation is introduced in 'Simulation'. The in- 
vitro experimental results are shown in 'In-vitro experiments'. The conclusions of this 
work are drawn in 'Conclusion'. 

Theory and method 

Photoacoustic theory 

In this paper, the two-dimensional PAI is concerned in the simulations and experiments. 
In 2D PAI, a laser pulse is used to illuminate the imaging tissues from the top. Due to the 
photoacoustic, the illumination creates an initial acoustical pressure field. The initial 
acoustical pressure field propagates as ultrasound waves, which can be detected by ultra- 
sound transducers. Based on the physical principle of the photoacoustic effect, assuming 
that the illumination is spatially uniform, the relationship between the acoustical pressure 
measurements and the initial pressure rise distribution can be derived as : 



where pyr ^tj is the acoustic pressure measurements at the position r and the time t, 
c is the sound speed, Cp is the specific heat, \i is the isobaric expansion coefficient, I{t) 
is the temporal profile of the laser pulse and ^i^^ is the initial pressure rise distribu- 
tion. In our study and many photoacoustic tomography studies, we employ a laser 
pulse with a very short duration. Its duration is nano seconds. So here we made an ap- 
proximation to treat the I{t) as a Dirac-delta function. 




(1) 
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In order to recover initial data for the wave equation, several inversion formulas have 
been established to solve this as a filtered back-projection problem [12,40]. By using 
the Green's function [12] to solve equation (1), the acoustic pressure measurements 
can be deduced as: 



r t 

o' J 4:71 Cpdt 



§LJ "-^^'7 (2) 



where r is the position of the ultrasound transducer. 

In PAI experiments, an ultrasound transducer is used to receive the acoustic pressure 
measurements at different positions, and the image reconstruction is regarded as an in- 
verse problem to obtain the initial pressure rise distribution. In the iteration of the 
image reconstruction, a projection matrix A is typically established to connect the 
acoustic pressure measurements with the reconstructed image. The measurements can 
be calculated based on the reconstructed image, and then the reconstructed image can 
be repeatedly corrected by minimizing the difference between the calculated measure- 
ments and the real ones. In this way, the optimization method can be used for collabor- 
ation and then the iteration reconstruction algorithm can be developed. 

Compressed sensing for PAI 

If the sampling data is insufficient, the projection matrix A is ill-conditioned. Thus, the 
matrix A does not have an exact inversion. As a result, it leads to streaking artifacts in 
the reconstructed image. This problem can be treated by incorporating the compressed 
sensing theory into PAI. 
We define a new variable /as: 

Then the equation (2) can be converted as follows: 

f{l^')=§i.-r,M'^yr (4) 

In practical imaging, the reconstructed image and the measurements are processed 
discretely, and the image is reshaped into vectors for convenience. If the size of the re- 
constructed image is X pixels x Y pixels, then the total pixel number of the re- 



constructed image "(^r j is N {N = XY). After vectorization, the reconstructed image 

M^7^ becomes a vector u with the length of N. If the total number of the detection 

points is Q, the length of measurement in each detection point is M, the equation (4) 
can be expressed as: 

/.=V.u i=l,2,--,Q (5) 

where / is the integration of the u along the arc that is centered in ith detection 
point and with a radius of ct, A, is the projection matrix of the ith detection point, 
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T is the transpose operation of a matrix. The calculation of the projection matrix 
is as follows: 



(a) Calculate an matrix Aiij) as: 



,0 (1<;<M) (6) 



where d = ^ [m-niif' + [n-riif' , {m, n) is the position of the /th point in the recon- 
structed image,(w„ «,) is the position of the iXh detection point,dc« is the actual length 
between the two pixels in the reconstructed image.di is the discretized time step and M 
is the total sampling points at one detection point. 



max< 1 



d ■ Ax 
c-dt 



(b) Vectorize the matrix Aiij) as the y'th column vector in projection matrix /I,. 

(c) Repeat the calculation M times to get the projection matrix A,. 

(d) Repeat step(a) to step{c) Q times to get the projection matrix in the different sampling 
positions {Ai,A2,...Aq). Then write the projection matrixes in the forms as follows: 



A= ^ 

w) 

The equation (6) can be expressed as: 

f = A-u 



(7) 



(8) 



where the sizes of / ^ and u are MQ pixels x 1 pixel, MQ pixels x N pixels and jV 
pixels X 1 pixel respectively. 

To reconstruct the photoacoustic image from incomplete measurements by using the 
compressed sensing theory, we can solve an optimization problem as follows: 



iin„||¥Tu|^ + ^|^u-/[2'| 



(9) 



where ^ is a sparse transform matrix, | |i and | I2 are the Li-norm and L2-norm re- 
spectively. By projecting the image onto an appropriate basis set, we can get a sparse 
representation of the original image. In this domain, most coefficients of the image are 
small, and a few large coefficients capture most information of the signal. In this way, 
we can recover a much more accurate image from those undersampled measurements. 

In practical applications of PAI, the reconstructed images often show piecewise con- 
tinuous behavior. The images like this always have small total variation (TV) values, 
which is defined as follows: 



P N 
TV{u) = / |r>u| = Ml i = 1. 2...iV 



(10) 



where D, is a matrix with the size of 2 pixels x N pixels that has two nonzero entries 
in each row to calculate the finite difference of u at the Jth pixel. D is a matrix with the 
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size of 2N pixels x N pixels, and D = {D^;D^),D^ and are the horizontal and vertical 
global finite difference matrixes respectively. 

It is reported that the TV based reconstruction algorithm can recover the image ac- 
curately from sparse sampling data [34]. Using TV values to reconstruct the image can 
be expressed mathematically as: 

min„|ri/(u)+^|Au-/|/| (11) 

However, the TV minimization still has some limitations that impair its performance. 
The optimization of the TV value encourages the recovery of images with sparse gradi- 
ents, thus resulting in the paintinglike staircase artifacts in the reconstructed images. 

Recently, some research find out that the nonconvex optimization can reconstruct an 
accurate image with fewer measurements by replacing the Li-norm with the L^-norm 
{Q <p < 1). Aiming to enhance the reconstruction quality and overcome the problem of 
TV based algorithm, we joint the L^-norm with TV values to establish a new opti- 
mization which can be defined by: 

min„ja.rnu)+;g|^T^^|/ + ^|^u-/|2^| 0<p<l (12) 

where a and |3 are parameters corresponding to the weights of the TV value and 
i^-norm value , | \p is the Lp -norm in this optimization problem respectively. 

Therefore, we can obtain the reconstructed image by solving this new optimization 
problem in equation (12). 

PAI reconstruction algorithm 

In this part, we solve the optimization problem in equation (12) to establish a novel 
photoacoustic image reconstruction algorithm by using the total variation and noncon- 
vex optimization. 

We define the finite difference approximations to partial derivatives of u at the it\\ 
pixel along the coordinate as variable o), = /),u, the ith pixel's sparse coefficient as vari- 
able z, = Wju, where is the sparse transform matrix of the ith pixel. The equation 
(12) can be deduced as: 

min,,,,„ |a^|c^,|2 + ^^|z,.|^ + ^|A.u-f|^| i=L2...iV (13) 

where p is the parameter corresponding to the weight of the constraint condition in 
this optimization problem. 
We form the augemented Lagrangian defined by 

(14) 

where &f is the TV step parameter in /rth iteration, cf is the L^-norm step parameter 
in kxh iteration, u'^ is the vectorized image reconstruction in ^h iteration. 
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This problem can be solved by 

' {(o'+\z'+\u^+') = min,,,,ui(6.,z,u;M,/), 

< b';+' = b^-{co^+'-Diu'<+'), (15) 
,cf+i =cf-(zf+i-'Fy), 

where 0,"*' is the finite difference approximations to partial derivatives of u in 
{k + l)th iteration, z'^^^ is the sparse coefficient in {k + l)th iteration, u**^ is the vector- 
ized image reconstruction in (k+ l)th iteration, zf^^ is the sparse coefficient of the ith 
pixel in {k+ l)th iteration, cof*^ is the finite difference approximations to partial de- 
rivatives of u at the ith pixel along the coordinate in {k+ l)th iteration; ^f^^ is the TV 
step parameter in {k+ l)th iteration, cf^^ is the L^-norm step parameter in (k+ l)th 
iteration. 

By using the standard augmented Lagrangian method, the optimization problem in 
(15) can be deduced as 

(16) 

where co'^ is the finite difference approximations to partial derivatives of u in 
^h iteration, z*^ is the sparse coefficient in /cth iteration, u^ is the vectorized 
image reconstruction in Ath iteration; 6^^ is the Barzilai-Borwein step parameter 
in Ath iteration. 

After using the Barzilai-Borwein method to determine the step size 6, the 
optimization problem in equation (13) can be transformed into three sub-problem 
as follows: 

' = min„,| + ^ - Au' - b1\l+^\co, - 

= min,{|..|;+^ |z, - ^Ju'< - cf |;+| |z, -41^}, 

< u'<+i= minuj ap\Du-co^+''\l+^p\'F'^u-z^'+'^\l+5k\u-{u''-S-^^A^{Au''-f))\lY 
^f+i = &f-(4+i-Au*+i), 

^ Sk+i = \A{n''+'-n%l (\co^+'-jf^+\z^+' - z'f^+\n^+' - u'f^) 

(17) 

where cof is the finite difference approximations to partial derivatives of u at the ith 
pixel along the coordinate in l<th iteration respectively, zf is the sparse coefficient 
of the ith pixel in kth iteration respectively; 6,^+1 is the Barzilai-Borwein step par- 
ameter in {k + l)th iteration. 

We use the soft shrinkage operator to obtain the solution to w-subproblem in equation 
(17), the operation is as follows: 
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< 



coM^ = max 



fli = Diu'' + bi 

k 



<2l + «2 



||l/(fll+fl2)| 



<a;2 = 
ti=p 

, i2 = 



{i = ia-N) 



(18) 



where Ai, iS!2« h and i2 are the variables used for a succinct expression. 
As for the z-subproblem in equation (17), we use the soft /7-shrinkage operator to 
solve it. The operator is defined by: 



.k+l 



ait 3 + a^ti 



Si + an 



a3 + Ui 



a-3 + fl4 



i-p 



l/(«3 + Ui) 
|l/(«3 + Ui) 



t-i= p 

ti = Sk/fi 



(i=l,2...N) 



(19) 



where 123, 124, and t^^ are the variables used for a succinct expression. 
The M-subproblem in equation (17) is a typical least squares problem. The solution 
can be easily obtained by: 



(20) 



k+i _ T/^[«/''DTQ^+i+^p<FTz'^+i + 5i:u'^-A-i(Au*-f) 

where F is the Fourier transform matrix. 

As a result, the TV- Lp algorithm is summarized as follows: 



(1) Initialization: input / a, |3, e ,p and p. Set the reconstructed image u° = 0, ft = c = 0, 
5o = l,k=0. 

(2) Apply equation (18) and (19) to update the value of q and z. 

(3) Apply equation (20) to update the value of u. 

(4) Apply equation (17) to update the value of b, c and 6. 

(5) If the exiting condition is met, end the iterations and output the result. Otherwise 
repeat the step from (2) to (4). The exiting condition is as follows: 



,k-l 



< e 



(21) 



Simulation 

To verify the effectiveness of the proposed TV- Lp algorithm on PAI reconstructions, 
the simulations are designed. All the simulations are performed in Matlab v7.14 on 
a PC with a 3.07 GHz Intel Xeon processor (only 1 core is used in computation) 
and 32 GB memory. The sparsisfying operator W is set to Haar wavelet transform 
using Rice wavelet toolbox. The sound speed is set to be consistent in the simula- 
tion as 1500 m/s. 
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Sparse-view reconstruction 

In the simulation, we choose the Shepp-Logan phantom to be the initial pressure rise 
distribution. The forward simulation and inverse reconstruction are all performed in 
2D. The phantom is shown in Figure 1. The measurements from the phantom are gen- 
erated by using equation (2). The size of the phantom is 89.6 mm x 89.6 mm, the radius 
of the scanning circle is 42 mm and the size of the reconstructed image is 128 pixels x 
128 pixels. During the simulation, the scanning circle covers 360°around the imaging 
phantom. Four different measurements are collected. The scanning step of tomo- 
graphic angels is set to 2.25°, 4°, 12° and 20° respectively. So the sampling points are 
160 views, 90 views, 30 views and 18 views correspondingly. 

The parameters a, p, e and p are set to be 1 x 10"^, 1 x 10"^, 1 x 10"^ and 1 respec- 
tively. The influence of these parameters will be discussed later. And the parameter p 
are set to be two different values as 0.5 and 0.8. 

We choose the FBP [16], the Li-norm [27] and the TV-GD [34] algorithms to be the 
comparison besides our proposed TV- Lp algorithm. The simulation results by using 
these different algorithms are shown in Figure 2. It's worthwhile to note that the weight 
used in the TV-GD algorithm is an adaptive parameter, as same as it is reported in 
[27]. The negative values in FBP reconstructed image are set to be zero. 

It is shown in the first column of Figure 2 that, all three iterative algorithms have 
comparable reconstruction results when the sampling data is sufficient. Moreover, it is 
shown in Figure 2(a) that the contrast of FBP reconstructed image is not as high as 
the other three. But its resolution is comparable with the others visually. When the 
number of sampling points reduces, the qualities of the reconstructed images are 
strongly affected in the FBP reconstruction. When the sampling point gets sparse in 
the FBP reconstruction, the arc-like artifacts appears due to the back-projection arcs 




Figure 1 The shepp-logan phantom. 
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Figure 2 The shepp-logan phantom reconstruction results by different algorithms. The first to fifth 

rows refer to the FBP (a-d), LI -norm (e-h), TV-GD (i-l), TV-Lp (p = 0.5) (m-p) and TV-Lp (p = 0.8) (q-t) 
reconstructed images respectively and the first to fourth columns refer to the results from 160-view, 
90-view, 30-view and 18-view respectively. 



cannot be canceled out with each other. The iterative algorithms can provide better 
qualities of the reconstructed images than the FBP method in sparse-view reconstruc- 
tions. Among them, the Lj-norm method struggles to depress the noise. Meanwhile, 
the TV-GD algorithm and the TV-L^ algorithm provides high-resolution images and 
have no visually distinguishable decline in qualities of the reconstructed images as the 
number of sampling points decreases. 

As for the extreme sparse sampling points situation (18-view and 30-view), the image 
reconstructed by the FBP algorithm, shown in Figure 2(d), has extremely severe artifacts. 
The Li-norm reconstruction and the TV-GD algorithm have a decline in image qualities. 
The noise in the reconstructed images by the Li-norm reconstruction, as shown in 
Figure 2(j), cannot be depressed effectively. As for the TV-GD algorithm, the reconstruc- 
tion produces piecewise artifacts which also make the qualities of the reconstructed 
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images decrease. In the TV-Lp image, the noise is depressed more effectively than the 
above three algorithm. The quality of the reconstructed image is not affected substantially 
by the insufficient sampling data. 

We calculate the peak signal-to-noise ratios (PSNR) of the reconstructed images with 
the original phantom as a gold standard to provide a numeric quantification of the re- 
sults. The bigger the value of PSNR is, the better quality of the image is. The PSNR is 
defined as: 



PSNR= lOlogi, 



XY -MAX 



X.Y 

\('V)=i / 



(22) 



where t{i,J) means the gray-value of the original image, MAX the maximum possible 
pixel value of the image which in our simulation is 1. 

We calculate the value of the PSNR of all images in Figure 2. The quantitative results 
are shown in Table 1. 

From Table 1, it is shown that the PSNR of the FBP algorithm is always in a very low 
level due to its unsuitability for sparse-view sampling condition. As for those three 
compressed sensing based algorithms, the PSNR value of images reconstructed by the 
TY-Lp algorithm are the highest. The L^-norm optimization constraint can provide the 
better performance in the extremely sparse sampling. With this improvement, the 
TY-Lp algorithm is more accurate than other algorithms in the sparse-view sampling 
condition shown in the quantitative results. Between the two different value of p, the 
parameter p that is set to be 0.5 has a slightly advantage against the other one. Also it 
is revealed from Table. 1 that the 90-view shows higher PSNR than 160-view. When the 
sampling points are sufficient, it is possible that the fewer-view projection can produce 
better reconstruction results. But their PSNR is very close with same algorithm. It is 
fair to say that the results are on the same level of image quality. 

From Figure 2, it is shown that the TV-GD image and the TY-Lp image are very close 
in the image quality. Here we choose the FORBILD phantom [41], a more complicated 
and more challenging phantom, to further compare the proposed algorithm with the 
TV-GD algorithms. The phantom is shown in Figure 3. The scanning step of tomo- 
graphic angels is set to 2.25°, 4°, 6° and 12°. So the sampling points are 160 views, 90 
views, 60 views and 30 views correspondingly. The other numerical implementation 
conditions remain the same with the shep-logan simulation. The simulation results by 
using the TV-GD algorithms and the proposed TV- Lp algorithm are shown in Figure 4. 
It is shown in Figure 4 that, when the sampling data is sufficient, both algorithms can 



Table 1 PSNRs (dB) of reconstructed images of shep-logan phantom 




1 60-view 


90-view 


30-view 


1 8-view 


FBP 


20.35 


18.36 


15.68 


13.14 


/_,-norm 


33.83 


34.98 


34.21 


31.19 


TV-GD 


38.01 


38.23 


36.68 


34.68 


IV-Lpip = 0.8) 


38.45 


39.05 


36.91 


36.72 


lV-Lpip = 05) 


38.85 


39.27 


37.01 


36.81 
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reconstruct the accurate image. When the sampling angles get sparse during the simu- 
lation, it is seen that the TV-GD reconstruction results in paintinglike staircase artifacts 
in the smooth regions. Also it fails to give the accurate image in the low contrast re- 
gions in the top and left of the phantom. The proposed algorithm provides reasonably 
good reconstructions in these regions. The PSNR of the reconstructed images are 
shown in Table 2. From this table, we observe that the TV-L^ algorithm provides the 
better PSNR for all the cases. In the case of more complicated phantom, the TV-L^ al- 
gorithm shows significant improvement than the TV-GD algorithm. 

Also, we include a line-plots image of the reconstruction result by the TV-GD algo- 
rithm and the TV-Lp{p = 0.8) algorithm from 30-view data. The location of the pixel 
profile in the image is displayed in Figure 5(a). The comparisons of pixel profiles are 
displayed in the Figure 5(b). 

In Figure 5(b), the solid line and the dotted line represent the pixel profiles of the 
TY-Lp and the TV-GD image respectively. It is shown in Figure 5(b) that the TV-Z.^ 
can reconstruct the image more precisely than the TV-GD one. The edges from the 
TY-Lp are sharper than that from the TV-GD. The pixel number from 90 to 100 is the 
high resolution area, the TY-Lp image shows the high-speed change of the pixel value 
while the TV-GD fails to do so. In the continuous area, the TY-Lp image is smoother. 

We continuously decrease the number of the detect points try to find the limit dens- 
ity of the sampling points. During the simulation, we set the criterion of acceptability 
that the PSNR of the reconstructed image reaches 30 dB. It is found out that the total 
number of the sampling points is able to be reduced to 15 for TY-Lp algorithm in the 
reconstruction of shep-logan phantom and to 18 for the forbild phantom. 

In this part, the TY-Lp algorithm is proved to be more accurate and stable than the 
other algorithms for PAI image reconstruction in the sparse sampling condition. 

Convergence and calculation 

In this part, we discuss the theoretical calculation complexity and study the conver- 
gence of the proposed algorithm. As mentioned above in 'Theory and method) in step 
(2) the update of co and z is using the soft shrinkages and the computational costs are 
both 0(7V). The update of z also includes a wavelet transforms which computational 
costs are 0(MogA/). In step (3), the update of u involves two fast Fourier transforms 
which computational costs are 0(MogA/) and two operations of A with the computa- 
tional costs of 0{NMQ). The update of the parameters b and c in step (4) are all simple 
calculations with the computational costs of 0{N). As for the paremeter 6, although it 
involves an operation of A, it can use the result computed in the step (3). So its compu- 
tational cost is also 0{N). 

In a nutshell, the calculation complexity of the proposed algorithm in one iteration 
is 5O(A0 + 40(Mog7V) + 20(NMQ). The first two terms is much smaller than the last 
term in the practical use of photoacoustic imaging and most iterative algorithm is with 



Table 2 PSNRs (dB) of reconstructed images of FORBILD phantom 





1 60-view 


90-view 


60-view 
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TV-GD 
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26.63 


IV-Lpip = 0.8) 
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Figure 5 The gray-scale profiles, (a) The location of the pixel profile in the image. (b)The gray-scale 
profiles of the reconstructed Images from 30-view data. 



the operation. In each iteration, we just use the projection matrix twice. So the pro- 
posed algorithm has a cheap per iteration computation. 

The TV-GD algorithm is reported as an efficient and stable iterative algorithm in 
photoacoustic imaging. In the 'Sparse-view reconstruction! its reconstruction result is 
closest to the proposed algorithm. So here we select it to be a comparison with the 
TY-Lp algorithm. We calculate the time cost of those two algorithms in a simulation. 
The simulation condition is same as in 'Sparse-view reconstruction'. But the iteration 
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ends when the PSNR values reach 30 dB. The result is shown in Table 3. From Table 3, 
it is shown that the proposed algorithm is faster than TV-GD algorithm in the compu- 
tational time. Based on this result, it could be inferred that the TY-Lp algorithm is a 
more efficient image reconstruction algorithm comparing to the TV-GD algorithm. 
The value of p has also some influence on the time cost. The smaller the p is, the more 
iteration times are needed to reach the reconstruction result. 

Thanks to the use of Barzilai-Borwein step size selection method, the convergence 
speed can also be significantly improved. For the quantitative analysis, we use a param- 
eter that represents the distance between the reconstructed image and the original 
phantom image. The parameter d is defined as: 

d= 

X Y 

where u is the reconstructed image and t is the original image. The size of the image 
is Xx Y. The smaller the parameter d is, the closer the reconstructed image is with the 
original phantom. In the TV-L^ algorithm, there is a small rate of chance that the 
optimization will lead to the wrong solution due to its non-convex nature. So we use 
the original image to calculate the parameter d to show the image quality and use the 
parameter d as a reference. We want to show the improvement of the image quality in 
every iteration step. 

The simulation condition is set to be the same as in 'Sparse-view reconstruction'. 
The sampling view is 60. The parameter p is set to be 0.8. The defined distance d is cal- 
culated after each iteration step. If the distance is smaller than 0.05, the iteration will 
stop. The simulation result is shown in Figure 6. The x-axis is the value of distance and 
the j-axis is the iteration times. The line '-' refers to the TV-GD algorithm and the line 
'*-' represent the TV-i^ (p = 0.8)algorithm. The result is shown in Figure 6. The images 
reconstructed by TV-L^ algorithm in each iteration have smaller value of d than the 
TV-GD ones and the TV-Z-^ iteration only takes 9 times as the distance is met the 
request. 

For discussions noted above, it can be surmised that the convergence of TV-L^ algo- 
rithm is faster and the TV-Lp has a cheaper computational cost. 

Robustness to the noise 

In the practical applications of the photoacoustic tomography, the measurements are 
usually polluted by those white measuring noises from the ultrasound transducer and 
the system electronics. Hence, it is very important for an algorithm to maintain stable 
performance under the noises polluted circumstances. To analyze the robustness of the 



Table 3 Time cost(second) of reconstructed images 





1 60- view 


90-vlew 


30-view 


1 8-view 


TV-GD 


41.79 


3719 


24.63 


14.61 


lV-Lp{p = 0.8) 


18.45 


12.05 


5.29 


441 


lV-Lp{p = O.S) 


20.76 


14.01 


8.69 


6.16 
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Figure 6 The distance between the reconstructed images and the original phantom image versus 
the iteration number. 



TY-Lp algorithm, we choose 30-view simulated photoacoustic signals that we used in 
'Sparse-view reconstruction'. The signals are added with white noises of different noise 
power levels. We use TY-Lp algorithm with two different settings of parameter p (p = 0.5 
and p = 0.8) and TV-GD algorithm to reconstruct images from these white noise polluted 
measurements. 

The reconstruction results are shown in Figure 7. In the first row to the last row of 
the Figure 7, the signal to noise ratio (SNR) of the polluted measurements is 10 dB, 
5 dB, 3 dB and 0 dB, respectively. As shown in the image, when the power level of the 
noises is not very strong (10 dB and 5 dB), the reconstructed images by using the noisy 
measurements have basically no obviously difference with the ones reconstructed with 
the noiseless signals. As the noise becomes stronger, the quality of the reconstructed 
images decreases. 

We also plot the profiles of a pixel line in order to show the detail qualities of recon- 
structed images clearly. In Figure 7, the dotted line and the solid line are the pixel pro- 
files of the reconstructed image and the original image, respectively. From the line 
plots, it is revealed that the proposed algorithm has better performance in the edge 
preservation and more accurate in the smooth area. We calculated the PSNR of the re- 
constructed image. The result is shown in Table 4. Our algorithm outperformed the 
TV-GD algorithm in any noise power level. Giving the credit to the optimal conditions, 
the reconstructed image is intended to be continuous and sharp. During the iteration, 
the photoacoustic signals get enhanced and the noise is suppressed. 

As we can see from the table that the TV-L^ algorithm reconstructed images have a 
slightly better performance than the TV-GD algorithm when it comes to the image 
qualities. When the noise is extremely strong (OdB), the TV-L^ algorithm has a huge 
advantage than the TV-GD one in image quality. As for the different settings of para- 
meter p in the TV-i^ algorithm, there is no major difference that can be observed be- 
tween the two images. The PSNRs of the p = 0.5 setting is about 0.3 dB bigger than the 
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Figure 7 The reconstruction results from noisy signal. The first to the third columns show images 
reconstructed by the 7V-GD (a-d) and 7V-Lp algorithm with p = 0.8 (e-h) and p = 0.5 (i-l), respectively, from 
signals with SNR of 10 dB, 5 dB, 3 dB and 0 dB. 



p = 0.8 setting in the first three noise power level. But when the noise getting extremely 
strong (0 dB), the p = 0.8 setting is 0.1 dB bigger than the p = 0.5 setting in PSNR value. 

From this part of simulation we can conclude that our TV-Z,^ algorithm is robust to 
noise and has a better performance than the TV-GD algorithm in the noisy measure- 
ment circumstance. 



Parameter investigation 

As the original optimization problem of the image reconstruction is described in Eq. 
(14), the TY-Lp algorithm contains some parameters that are tunable, which are a, |3, e, 
p and p. In those parameters, the choice of p does not affect the performance of the 



Table 4 PSNRs (dB) of the noise simulation image 





10 dB 


5 dB 


3 dB 


0 dB 


TV-GD 


32.24 
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22.44 
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lV-Lp{p = 0.8) 


35.14 
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25.21 


TV-Lp(p = 0.5) 
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28.10 


25.06 
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TV-Lp algorithm theoretically. The result of the simulation also shows that the image 
quality is not sensitive to the parameter p for a large range. Here we set the parameter 
p to a steady value 1. The e is the exiting condition parameter. It can be easily deduced 
that smaller e will leads to slightly more accurate reconstructed image at the cost of 
more iteration times. In this part we focus on analyzing the parameter settings of p, a 
and p. 

Parameter setting of p 

In the TY-Lp algorithm, we replace the Li norm with the Lp norm (0 <p < 1). It is re- 
ported in Ref [35] that theoretically fewer measurements are required for accurate re- 
construction in the Lp norm situation. But it also leads to failure in solving the 
optimization problem. It's kind of a dilemma for the setting of p. So here we take differ- 
ent values of p to see its influence to the image reconstruction. The parameter a and (3 
are both set to be 1 x 10"^. 

We choose the 90-view and 18-view simulated photoacoustic signals that we used in 
'Sparse-view reconstruction'. We set the p value as 0.3, 0.5, 0.8 and 1. We calculate the 
reconstructed images' PSNR value. It is shown in Table 5. When the p is set as 0.5, it 
has advantage in the quality of the reconstructed image. However, when the value of p 
continues to reduce to 0.3, there is no obvious improvement in image quality. But in 
the same time, the smaller the p is, the higher probability of the solving failure is dur- 
ing the simulation. The reduction of p leads to the increasing of the iteration times in 
our simulation. So taking these two factors into account, we set the p as 0.8 so that it 
can provide a great reconstruction performance and stability with a fast convergence. 

Parameter settings of a and p 

As we describe above in Eq. (14), the parameter a and p are parameters corresponding 
to the weights of TV value and Lp norm value in this optimization problem respect- 
ively. We use these two parameters to balance the terms of the objective function. With 
different kinds of the objective image, the settings of those two parameters are differ- 
ent. Here we select three different images as the given optical energy deposition to test 
the universality of our algorithm and present further investigation of the parameter set- 
tings. We select a phantom that stand for the vessels and a phantom of dots with dif- 
ferent energy degree in the simulation. We also choose a real brain MRI as the original 
optical energy deposition to demonstrate the performance of the proposed algorithm in 
reconstructing extremely detailed and complex structured imaging object. Here the 
TV-GD algorithm is used as a comparison. There are four groups of parameter settings, 
which are (a = 1 x 10"^, p = 5 x 10"^) , (a = 1 x 10"^ p = 1 x 10"^), (a = 5 x 10"^, P = 5 x 
10 and (a = 5 X 10"^, p = 5 x 10"^). The reconstructed images are shown in Figure 8. 
From first row of Figure 6, we can see in the reconstruction of gradient sparse phantom. 



Table 5 PSNRs (dB) of reconstructed images with different p 
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Figure 8 The different image objectives reconstructed by TV-GD and TV-Lp algorithms. The first to sixth 
rows refer to the original phantom (a, g, m), the TV-GD reconstructed images (b, h, n), the TV-Lp reconstructed 
images with parameter settings as (q = 1 x 1 0~^, p = 5 x 1 0"^) (c, 1, o), (a = 1 x 1 0^^, (5=1x1 0"^) (d, j, p), 
(a = 5 X 10"^ P = 5 X 10"^) (e, k, q) and (a = 5 x 10"^ p = 5 x 10"^) (f, I, r), respectively. 



The TV based algorithm has great performance when the image demonstrate piecewise 
continuous behavior. All reconstructions are accurate and the background noise is sup- 
pressed well. When it comes to the images with the vessel phantom (Figure 8 (g)-(l)), 
those original optical energy depositions are a little bit more complex than the dots. The 
reconstruction results show that the image reconstructed by TV-i^ algorithm is better 
than the TV-GD ones. As in Figure 8(h), TV-GD image has some noises in the back- 
ground and the edge of the vessel is blurred. While the TV-Lp images with different par- 
ameter settings both have high-resolution results. As for the real MRI image, it has very 
detailed information. As expected, both two groups of parameter setting a = (3 have the 
most accurate result among them. The increasing weighting of L^-norm condition can 
provide more detail information and prevent the reconstructed image emerging plantlike 
artifacts. The details such as edges and fine structures are well preserved in both recon- 
structions. The reconstruction results show that the TV-GD reconstructed image has se- 
verely paintinglike staircase artifacts with some loss in fine details. From our observation, 
TV-Lp algorithm with the parameter setting a = (3 preserves the fine features better than 
the TV-GD one. a and |3 are the regularization parameters determining the trade-off 
between the data consistency and the sparsity. It is revealed from the above simulation 
that the parameter setting a = p is a better strategy. In this parameter setting, the TV-Lp 
algorithm provides a 3 dB improvement in the PSNR over the TV-GD algorithm based on 
our calculation. 

Limited-view and irregular-view simulation 

In the real application of PAI, due to the restrains of the shape or the size of the 
imaging object, a full angular scanning sometimes is hard to achieve. We evaluate 
the performance of the TV-Lp method in limited-view case, line-view case and un- 
equal-view case. 
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The simulation setup and the reconstructed image is shown in Figure 9. In the 
Umited-view simulation (Figure 9(a)), the scanning angular range is set to 150° and the 
angular step is 3°, so 50-view photoacoustic signals are obtained. In the line-view simu- 
lation (Figure 9(c)), the transducer array with 60 transducers is placed in the right side 
of the imaging object and the interval between two transducer elements is 1.49 mm. It 




(e) (0 

Figure 9 Limited-view, line-view and un-equal angel step scanning reconstruction results, (a) refers 
to the location of transducers in the limited-view simulation, (b) refers to the reconstructed image from the 
limited-view data, (c) refers to the location of transducers in the line-view simulation, (d) refers to the 
reconstructed image from the line-view data, (e) refers to the location of transducers in the un-equal angel 
step scanning, (f) refers to the reconstructed image from the un-equal angel step scanning data. 
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is revealed in Figure 9(b) and Figure 9(d) that the quahty of the TV-Lp reconstruction 
is not much affected by the limitation of the sampling angle. Because the sampling 
angle is limited, the information definite is partly missed, yet the TY-Lp method can 
still provide a satisfying reconstruction. In un-equal angel step scanning, we randomly 
choose 30 sampling points from a 60-view projection and use these 30-view un-equal 
angle step data for image reconstruction. The result is shown in Figure 9(f). As we can 
see from the image, the reconstruction result can still maintain a very high quality. 



In-vitro experiments 

Experiment setup 

We carry out the experiments on in-vitro signals to demonstrate the proposed TY-Lp 
algorithm's performance in the practical application. 

The framework of the experiment platform is shown in Figure 10. In this platform, 
an Nd:YAG laser generator (Continum, Surelite I) is used to emit the laser pulse. The 
wavelength of the laser is 532 nm. A single laser pulse is generated at the frequency of 
10 Hz and last 6-7 ns. The incident laser pulse is emitted towards the top of the phan- 
tom through a concave lens with the diameter of 5 cm. The setup of the lens enlarges 
the illumination area and lead to the pulse energy reduction in the illumination area. 
The energy is about 6.47 mjcm"^, which is lower than the ANSI laser radiation safety 
standard (20 mjcm"^) [1]. Signal acquisition is done by a water-immersion ultrasound 
transducer (Panametric, V383-SU). The transducer is a linearly unfocused one at 3.5 
MHz (-6 dB bandwidth at 45%). A digital stepping motor (GCD-0301 M) is used to ro- 
tate the ultrasound transducer around the phantom placed in water. The scanning ra- 
dius is 38 mm. The received analog ultrasound signals are amplified by a pulse receiver 
(Panametric, 5900PR). An oscilloscope (Agilent, 54622D) with the sampling frequency 
of 16.67 MHz is set to transform the received signals into digital ones. Both the laser gen- 
erator and the digital motor are controlled by the computer through the serial interface. 
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Figure 10 The scheme of the experiment platform. 
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The transformed digital data is transported to the computer through the general purpose 
interface bus (GPIB). 

The imaged phantom we used in the experiment is made by gelatin cylinder. It is 
shown in Figure 11. There are two different phantoms. The radius of the phantom is 
25 mm. The left one is made by two rubber bars with 1 mm diameter that embedded 
as the optical absorbers. The right one utilizes leaf which pretends as vein and tissue as 
the optical absorbers. 

In the experiment, the transducer tends to measure the photoacoustic signal in-plane 
only, and the reconstruction is also in 2D. The cross-sectional image in any plane is 
mainly determined by the measured data in the same plane, and a set of circular meas- 
urement data on the same plane would be sufficient to reconstruct a good image. We 
use the deconvlution calculation before the reconstruction to eliminate the transducer's 
impulse response influence. 

Experiment result 

In the experiment, 90-view and 30-view data are collected for reconstruction. The im- 
ages are constructed by the FBP, the TV-GD and the TV-L^ algorithms, respectively. 
The reconstruction results are shown in Figure 12. The left column of the Figure 12 is 
reconstructed from 90-view data. When the sampling data is sufficient, all three algo- 
rithms are effective. With respect to the locations and sizes, the optical absorbers are 
all well reconstructed in the figure. While the FBP reconstructed image is not as clear 
as the images reconstructed by the iterative algorithms. When we reconstruct the 
image with a small of sampling angles (right column of Figure 11), the artifacts start to 
emerge in the FBP reconstructed image and the quality of the image is severely af- 
fected. But the TV-GD and the TV-Z.^ algorithms can still provide high-contrast images 
with less noise. In Figure 12(f), it is shown that the image reconstructed by the TV-Lj, 
algorithm outperforms other algorithms in image contrast and noise suppression. The 
structure of optical absorbers is clear and the noise in the background is well- 
suppressed. The sparse-view of sampling has barely any influence on the quality of 
the TV-Lp reconstructed image. 

In vitro imaging of a leaf vein is also performed to further demonstrate the advan- 
tages of the TY-Lp algorithm. The reconstruction result is shown in Figure 13. As the 




Figure 11 The photos of the imaging samples. 
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Structure of the phantom is more complex, the FBP are deeply influenced by the ar- 
tifact and fail to reconstruct the accurate image both under 90-view and 30-view sam- 
pling circumstance. It is shown in the figure, that the TV-GD and TV-Lj, algorithms 
can still reconstruct the image in a high contrast level. But when the data is insufficient, 
there is some noise emerging in the background. TV-L^ algorithms can suppress the 
noise better than the TV-GD one. The optical absorber in TV-L^ one is more distinct 
than that in the image by TV-GD algorithm. 



Quantitative comparisons 

We use the Z,i-norm algorithm to reconstruct the image of 180-view data from the leaf 
vein phantom. As the sampling view is efficient, the reconstructed image is used as a 
"standard" one. We calculate the histograms of the difference between the reconstructed 
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one and the "standard" one as shown in Figure 14. Figure 14 (a)-(c) are the difference his- 
tograms between the standard and the images reconstructed by the FBP, the TV-GD and 
the TY-Lp, respectively, with 30-view data. In Figure 14, two CS-based algorithms have a 
large number of pixels with small ranges of difference with the standard one, which sug- 
gests that these two algorithms can reconstruct the image more accurately. In the case of 
the TV-Lp algorithm, the major part of the pixel difference is in the range from 0 to 0.1. 
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Figure 14 (See legend on next page.) 
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f ^ 

(See figure on previous page.) 

Figure 14 The amplitudes histograms of the difference between reconstruction results and the 
"standard" one. We use the image reconstructed by Ll-norm with data from all 180 transducer elements 
as the standard, (a) Histograms image of the FBP algorithm, (b) Histograms image of the TV-GD algorithm, 
(c) Histograms image of the TV-Lp algorithm. 



From this experiment, the results demonstrate that the TV-L^ method can outperform 
the TV-GD one in the field of the image quality. 

From the experiment result noted above, it is safe to say that the TV-Z,^ algorithm 
would have better performance in sparse-view PAI than other algorithms. It could pro- 
vide stable and accurate reconstruction in both sufficient data sampling and sparse- 
view sampling situation. 

Conclusion 

Aiming to reduce the scanning time and enhance the imaging quality of the photo- 
acoustic image reconstruction, we proposed the TV-Z.^ algorithm that applies the total 
variation method and nonconvex optimization method to the PAI. The main idea of 
the algorithm is to apply L^-norm nonconvex optimization along with the total var- 
iation method. In the proposed algorithm, the Barzilai-Borwein step size selection 
method is adopted to provide faster convergence and smaller calculation. The effective- 
ness and universality of the algorithm is demonstrated through the numerical simula- 
tions. The numerical simulations show that the TV-Lp algorithm provides good imaging 
quality in sparse-view sampUng situation. The algorithm convergence, the robustness to 
noise and the tunable parameters are also discussed. The simulation result reveals that 
the TV-Lp algorithm is a stable image reconstruction method with fast convergence and 
small computational cost. The TW-Lp algorithm is further investigated through some ex- 
periments using gelatin-made phantom. Compared with the result of other popular image 
reconstruction method, the TY-Lp imaging algorithm has significant advantage on con- 
trast and noise suppression. From the discussion noted above, it could be concluded that 
the TY-Lp algorithm may be a practical algorithm for sparse-view photoacoustic imaging 
reconstruction. 
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